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■ In the framework of the LDA+U approximation we propose the direct way of calculation of crystal- 

' field excitation energy and apply it to La and Y titanates. The method developed can be useful 

for comparison with the results of spectroscopic measurements because it takes into account fast 
relaxations of electronic system. For titanates these relaxation processes reduce the value of crystal- 
field splitting by ~ 30% as compared with the difference of LDA one electron energies. However, 
the crystal-field excitation energy in these systems is still large enough to make an orbital liquid 
£SJ . formation rather unlikely and experimentally observed isotropic magnetism remains unexplained. 

PACS numbers: 71.15.-m, 71.20.-b, 71.30.+h 
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J I. INTRODUCTION 

The magnitude of crystal-field splitting (CFS) is known to be very important characteristic of transition metal 
C$ . compounds, necessary for understanding different physical phenomena such as magnetism and metal-insulator 
transitions^^ It is often used for detailed fitting of spectroscopic data and in various model calculations. An important 
question is: how can one calculate CFS in a reliable way? 

One can obtain character of CFS using the group theory analysis. It provides qualitative description of ionic levels 
structure in the presence of crystal-field of given symmetry and can even give some ratios between splittings using 
£j ' Wigncr-Eckart theorem^ However, for quantitative estimates the detailed structure of ionic potentials has to be taken 
into account. 

Generally, there are different contributions to the total value of CFS to be considered. A point-charge contribution 
can be calculated "by hand" using for instance direct Madelung potential summation (see e.g^). However, this 
^ r | . method has some drawbacks and neglects the overlap of the electron clouds of neighboring ions, which can lead to 
00 ' several important consequences i 

CN| , Another important contribution comes from the covalency. It can also be estimated using some model calculations 
(for instance taking Slatcr-Koster parameters), but without precise knowledge of the band structure one does not 
actually know the parameters of the models accurately enough. In order to take into consideration all these contribu- 
tions to the total value of CFS in a reliable way one should use ab initio methods based on density functional theory 
(DFT) calculations. Nevertheless, even in this case there are several ways to define and calculate CFS. Each of them 
has its own virtues, but also difficulties. 

The first possibility is a direct calculation of the center of gravity of corresponding bands. It works well when 
splitting is large enough. Thus it is most useful and rather accurate for t2 g — e g splitting (typically in 3d-oxides 
CFS between t2 g and e g states is ~2-2.5 eV). For example the behavior of spin-state transition temperature in 
cobaltites, which depends on this splitting has been recently estimated quite accurately^ The problems might arise 
O ■ in computation of smaller values such as the splitting inside t2 g or e g shells. The general problem is that CFS is a 
well-defined characteristic only for isolated ions which have energy levels but not bands. The bigger the corresponding 
band width (in comparison with CFS) the more questionable it becomes to treat the difference of centers of gravity 
IT1 [ as an estimate of CFS value. 

The second way to calculate CFS is the downfolding or projection procedure which could give a few orbital on-site 
Hamiltonian in the minimal basis set of functions ip 2 , V'JV^™ In Sec. IV the method of obtaining such few 
orbital Hamiltonian from a full orbital one is briefly described. 

Having this Hamiltonian one can diagonalizc it to obtain its eigenvalues ei, £2, £n and eigenvectors V^i, \&2> 
'Fjv which are in general linear combinations of the initial wave functions: 
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The eigenvalues can be considered as the energies of the orbitals. The differences between these eigenvalues define in 
this case CFS. 
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In spite of great generality of this method, it has some disadvantages. There is an ambiguity in the definition of 
unitary transformation matrix Uj^ for Wannier function construction (see Sec. IV) and in the choice of the basis 
set wave functions ip%, ipN in any downfolding or projection procedure. Using different basis sets few orbital 
non-interacting Hamiltonian can be constructed in different ways. All of them are equally good if the resulting 
Hamiltonian gives the same bands as the LDA full orbital one. It is essentially important for low symmetry systems, 
where it is not clear what kind of linear combinations of d— wave functions and in which local coordinate system 
(LCS) should be taken as a basis set for few orbital Hamiltonian construction. 

In addition there is one more serious and general problem. LDA calculations are widely known to give good 
description of the ground state characteristics. However, it often fails to describe the excited states. LDA eigenvalues 
should not be treated as the real excitation energies. This general shortcoming concerns all methods of calculating 
CFS which use LDA eigenvalues. 

On the other hand, care should be taken in what we actually mean by CFS. In direct study of CF excitations e.g. by 
optics one should take into account possible relaxation of both the electronic (fast process) and lattice (relatively slow 
relaxation) subsystems. Fast electronic relaxation definitely has to be taken into consideration, whereas the lattice 
relaxation can be often treated separately (direct optical absorptions usually occur at a frozen lattice, according to 
the Frank-Condon principle). 

In this paper we propose the direct method of calculation of CF excitation energy in the framework of the LDA+U 
approximation which allows for such relaxation of electron system and enables to avoid ambiguity in the choice of 
basis set (as in projection or downfolding procedure) as well. CF excitations are calculated as difference between the 
ground state energy and the total energy of an excited state in which the electron is "artificially" put to one of the 
higher-lying d— levels. These states - the ground state and the excited state with the electron constrained in the higher 
level - are both treated in the self-consisted LDA+U scheme. This allows us to obtain the values of CF excitation 
energy which, in particular, take into account fast relaxation of the electronic system. 

Below we develop this general method and test it on the example of Y and La titanatcs. We compare our results 
with the previously obtained ones&S The physics of these compounds is discussed on the basis of calculated values 
of CF excitation energy, as well as the other results of ab initio LDA and LDA+U calculations. 

II. PHYSICS OF TITANATES AND CRYSTAL-FIELD SPLITTING 

Unusual and very rich physics of titanates attracts much attention since the possibility of existence of an orbital 
liquid has been proposed for LaTi03ii2iii In the ground state both LaTiC-3 and YTiC>3 crystallize in pcrovskite 
structure. The lattice distortions due to different ionic radii of Y and La ions seems to lead to the cardinal changes 
in electronic and magnetic properties of these compounds. 

At low temperatures LaTiC>3 was found to be G-type antiferromagnet with Neel temperature for stoichiometric 
samples Tjy=146 while YTi03 is an isotropic ferromagnet with relatively low Tc ^30 Local magnetic moment 
on Ti 3+ ion in YTi0 3 was found to be 0.84 ll b ~ The ordered moment in LaTi0 3 amounts to 0.46 - 0.58 /x i/' 10 i 13 and 
strongly differs from 1 /j,b expected for d 1 configuration with quenched orbital moment. Another interesting feature 
is nearly isotropic magnon spectra with small spin gap observed both in LaTi03 and YTi03 despite of the intrinsic 
orthorhombic distortion>i2ii4 

In order to describe these puzzling properties the exciting idea of the orbital liquid has been proposed^ According 
to this theory the large degeneracy of ti g shell can lead to the orbitally-disordered ground state. The crucial point 
in this case is the value of energy required for electron excitation from one t2 g orbital to another: orbital liquid can 
exist only if this energy is zero or very small. 

However, as the real symmetry of LaTi03 and YTi03 is not cubic but orthorhombic, one should expect certain split- 
ting of ^2g-levels. First calculations of the CFS in LaTi03 was performed by R. Schmitz and E. Miiller-HartmannAi^ 
Using realistic crystal structure^ they carried out model calculations taking into account point-charge and covalcncy 
contributions and obtained that the lowest singlet is separated from two higher-lying almost degenerate levels by 
about 0.24 eV. Similar calculations were also performed by Mochizuki and Imada>i& They stressed the importance of 
the GdFe03-type distortion and found the value of CFS to be 0.77 /enLa eV, where enLa is an effective dielectric 
constant. This constant is hard to evaluate in model calculations in a reliable way because of local screening effects 
in a solid, which could be explicitly taken into account only in ab initio band structure calculations. 

The ab initio calculations were done by Pavarini et al& and Solovyev^ They used the diagonalization of few 
orbital on-site effective LDA Hamiltonian in real space to obtain CFS. Solovyev has performed an exact procedure of 
full-orbital Hamiltonian transformation to the small energy-dependent Hamiltonian. After that the energy 

was fixed in the center of the ti g bandit Pavarini et alA have used formalism of Nth-order muffin-tin orbitals (NMTOs) 
to define small Hamiltoniant£ Despite the similarity of the methods, the results arc qualitatively different. The CFS 
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between lowest and next t2 g level obtained in Ref. is rather small, ^30-50 meV for both compounds, which would 
not prevent an orbital liquid formation. However, the results of calculations presented in Ref. are incompatible 
with this scenario: obtained values of CFS are ~ 150-200 meV for both compounds, similar to the results of model 
calculations 



III. LDA: BAND STRUCTURE 



Discrepancy between the results of previous ab initio CFS calculations^ is probably caused by (i) choice of MT 
spheres radii in LMTO method^ or (ii) the details of different projection procedures. To resolve first problem we have 
carried out the conventional LDA calculations using LMTO method and verified its band structure by performing 
the full-potential calculations in the framework of linearized augmented plane waves (FP-LAPW) method realized on 
Wicn2k program code^ We chose FP-LAPW because it is known to be the most accurate method for band structure 
calculations. 

Crystallographic data for LaTi0 3 (T=8 K) and YTi0 3 (T=293 K) were taken from Ref.H The radii of the muffin- 
tin (MT) spheres for LMTO calculations were chosen to be RL a =3.28 a.u.. Rxi=2.52 a.u., and for both oxygens 
Ro=l-92 a.u. Remaining unit cell volume was filled by empty spheres (atomic spheres with zero nuclear charge) with 
different MT radii. Ionic radius of Y is known to be smaller than La one^ According to this fact the MT radius of 
Y was also taken to be smaller: Ry=3.01 a.u. 

The Ti(4s,4p,3eQ, 0(3s,2p.3<i) and Y(5s,5p,4d,4/), La(6s,6p,5<i) states were included in the basis set in our calcu- 
lations. Almost empty La-4/ states were treated as pseudo-core, since La-4/ states are localized and do not strongly 
contribute to the relevant electronic states. Due to rather large values of the on-site Coulomb interaction U (of the 
order of 10 eV) La-4f states will be located approximately at 5-10 eV above the Fermi levels 

The Brillouin-zone (BZ) integration in the course of the self-consistent iterations was performed over a mesh of 27 
k-points in the irreducible part of the BZ. Density of states (DOS) were calculated by the tetrahedron method with 
512 k-points in the whole BZ. 

The band structure obtained within the LDA approximation in the framework of LMTO (with the MT radii chosen 
above) and FP-LAPW methods for LaTiOs is presented in Fig. 2] Twelve Ti— tig (four Ti per unit cell) bands are 
placed in the vicinity of the Fermi level and have band width ~ 1.95 eV both in FP-LAPW and LMTO methods. 
The energy gap between the top of the ti g and the bottom of e g bands is estimated to be 0.25 eV in FP-LAPW and 
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0.18 cV in LMTO. Comparing this figure and Fig. 2 of Rcf. one can see that the present LMTO bands agree better 
in all high-symmetry points of BZ with FP-LAPW ones than those presented in Ref. 0. So being firmly convinced of 
the correctness of LMTO band structure one could carry out the CFS calculations. 



IV. LDA: WANNIER FUNCTION PROJECTION FOR CRYSTAL-FIELD SPLITTING CALCULATION 

One of the ways to calculate CFS is to use the minimal basis set Hamiltonian for the orbitals of interest in real 
space. To construct such Hamiltonian we used Wannier functions projection procedure. In this section we present 
brief description of the method. For more details, see Ref. 0. 

For an LDA Hamiltonian H one has a Hilbcrt space of eigenfunctions (Bloch states IV^k)) with the basis set \<j>^) 
defined by the particular method (for LMTO these are LMT-orbitalsji£ for LAPW - Augmented Plane Waves^I etc.). 
In this basis set the Hamiltonian operator is defined as: 

H = J2\M H »»(<t>»\, (2) 

l 1:1/ 

where greek indices are used for full-orbital matrices. 

If we consider a certain subset of the Hamiltonian eigenfunctions, for example Bloch states of partially filled bands 
l^rtk), then we can define a corresponding subspacc in the total Hilbcrt space. The Hamiltonian matrix is diagonal in 
the Bloch states basis, however, physically more appealing is a basis which would have a form of site-centered atomic 
orbitals. That is a set of Wannier functions |W^} defined as a Fourier transformation of certain linear combination 
of Bloch functions belonging to this subspace (see Eq. 10} bellow). They are labeled in real space according to the 
band n and the lattice vector of the unit cell T which they belong to. The Hamiltonian operator H WF defined in 
this basis set is 

H WF = J2 \WZ)H nn ,(T)(W%\. (3) 

nn'T 

The total Hilbert space can be divided into a direct sum of above introduced subspace (of partially filled Bloch 
states) and the subspace formed by all other states orthogonal to it. Those two subspaces are decoupled since they 
are the eigenfunctions corresponding to different eigenvalues. The Hamiltonian matrix in Wannier function basis is 
block-diagonal, so that the matrix elements between different subspaces in Hilbert space are zero. The block H nn i in 
© corresponding to the partially filled bands can be considered as a projection of the full Hamiltonian operator @ 
to the subspace defined by its Wannier functions. 

Localized Wannier functions |W i T ) were defined in Ref. |52 as Fourier transforms of the Bloch functions \ipiw) 

\ w ?) = ^E e " ikT i^>, ( 4 ) 

where N is the number of discrete k points in the first BZ. 

Wannier functions are not uniquely defined because for a certain set of bands any orthogonal linear combination 
of Bloch functions \ijjik) can be used in In general it means that the freedom of choice of Wannier functions 
corresponds to freedom of choice of a k-depended unitary transformation matrix 

3 

The most serious drawback of Wannier representation is that there is no rigorous way to define U^' . As one can 

see from (0} and (JSJ, Wannier functions can vary significantly in shape and range because variations in \ipik) or 
change relative phases and amplitudes of Bloch functions at different k and bands i. 

In order to avoid this disadvantage some additional restrictions on the properties of Wannier functions have been 
proposediSiSiSii Among others Marzari and Vandcrbildt— suggested the condition of maximum localization for Wan- 
nier functions. That gave the variational procedure to calculate . To get a good initial point the authors suggested 
to choose a set of trial localized orbitals \4> n ) and projecting them onto the Bloch functions It was found that 

this starting guess is usually quite goodi^ This fact led later to the simplified calculating scheme proposed inSi where 
variational procedure was abandoned and the result of aforementioned projection was considered as a final step. 
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TABLE I: The magnitudes of CFS in tig shell (in meV) obtained in the LDA approach by different authors are presented 
in first 3 columns. First value is an energy difference between the lowest energy level and the middle one; the second is the 
difference between the middle and the highest energy levels in t2 g shell. The results of the constrained LDA+U calculations, 
which take into account fast electron relaxations are shown in the last column. 
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In order to start projection procedure one needs to determine the set of trial orbitals \<j> n ) and the bands which will 
be used for the Wannier functions construction. The latter can be defined either by the bands numbers (from iVi to 
N2) or by the energy interval (E±, E%). 

Non-orthogonalizcd Wannier functions in real \W^) and reciprocal space jWrOc) are then the projection of the set 
of trial site-centered atomic-like orbitals \<p n ) on the Bloch functions l^k) of the chosen bands 



\WJ 



;L^e- ikT |M/ )lk ), (6) 



N k 



N 2 



\W nk ) = £ |V'ik)(V , ik|0n> = £ \i>ik}{lf>ik\<Pn)- 

i=JVi i(Bi<Ei(k)<_B 2 ) 

Note that the Wannier functions in reciprocal space jWrOc) do not coincide with the Bloch functions \ip n k) for 
multi-band case due to the summation over band index i in © . One can consider them as Bloch sums of Wannier 
functions analogous to the basis functions Bloch sums </>^(r), see (JSJ) below. 

The coefficients (tpiii\<t>n) in © define (after orthonormalization) the unitary transformation matrix in J5J. 
However, projection procedure defined in © can be considered as a more general formalism than unitary transfor- 
mation JSJ. 

The Bloch functions in LMTO basis take the form 



hfck) = (7) 



where /1 is the combined index of the qlm (q - atomic number in the unit cell, Im are orbital and magnetic quantum 
numbers), and <^(r) are Bloch sums of the basis orbitals </> M (r — T) 

#(r) = -^£e 2fcT ^(r-T), (8) 

and the coefficients have the property cj^ = (<^ 4 |i/ , ik)- 

If n in \4> n ) corresponds to the particular qlm combination (in other words \4> n ) is an orthogonal LMTO basis set 
orbital), then (ipik\4> n ) = c k * and hence 

N 2 

iww = J2 w*)<% ( 9 ) 

i=JVi 

l — Nl f-L f-L 



i=Ni 



k* 



6 



In order to orthonormalize the WF © an overlap matrix O nn ' (k) and its inverse square root S nn < (k) can be defined 

as 

N 2 

0„„,(k) = (W n *\W n ,*) = ]T C ™ C »^ ( 10 ) 

i=Ni 

s m ,Qc) = a:,y 2 (k), 

The orthogonality of Bloch states (V'nklV'n'k) = <W' was used. 

Orthonormalized Wannier functions in fc-space |W n k) can be obtained as 

N 2 

\W nk ) = E 5 -'(k)|W n 'k) = £ |Vik)c^ (11) 



i=JVi 



where 

-k* 

'■■■!■) i 



(V>ik|W„ k )=£Sw(k)c^, (12) 



iV 2 



= «|W«k) = £ i&cft. (13) 



Thus, matrix elements of the few orbital Hamiltonian H WF in the basis of WF in real space where both orbitals 
are in the same unit cell are given by the following expression 



HZ F (0) = (^ EE |^k}e i (k)(^k|)|W°} = 

^ k i=Ni ' 
N 2 

= ^EE^(%(%( k )' ( 14 ) 

k i=JVi 

Here £j(k) is an eigenvalue for a particular band. 

For both LaTiC>3 and YTiC-3 we are interested in CFS in the t2 g shell. For that 3x3 Hamiltonians were derived for 
t 2g bands placed in the vicinity of the Fermi level in the energy range (-0.8; 1.25) eV for LaTiC>3 and (-0.65; 1.55) eV 
for YTiC>3. Then we diagonalize them and calculate the difference between corresponding eigenvalues in order to 
determine the CFS. 

The obtained values of CFS together with results of Ref. 00 are presented in Tab. [I] The CFS calculated by 
Pavarini et alA has the same character and similar values as ours. According to these results, in LaTi03 the lowest 
energy level is widely separated from two other, which are almost degenerate. Similar result was obtained by Cwik et 
al. using a full Madclung-sum point charge model, where the first splitting was found to be 240 meVi 

The situation is a little bit different in YTiC-3. The magnitude of the first splitting (between lowest and middle 
energy levels) has the same order as in La titanate, while the second one (between middle and highest energy levels) 
increases. It is interesting to note that in calculations of Solovyev the same effect of the second splitting enhancement 
going from LaTiC>3 to YTi03 is observed^ Nevertheless, the character of CFS obtained by Solovye^ is quite different 
from the present one and from that calculated in Ref. 0. 



V. LDA+U: BAND STRUCTURE 



The application of LDA to solids provides important information about details of band structure. However, for 
transition metal compounds it usually leads to metallic type of the electronic structure due to the presence of the 
partially filled d— bands at the Fermi level, often in contrast with an experiment. This is also the case in the LDA 
approach presented above (Sec. IV). The LDA+U approximation has been developed to include in the calculation 
scheme orbital-dependent Hubbard-like correction U which acts differently on the occupied and unoccupied d— orbitals 
giving as a result correct description for localized states^ 
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FIG. 2: (color online). Total and partial DOS of LaTiC>3 (top panel) and YTi03 (bottom panel) calculated within the LDA+U 
approximation for AFM-G and FM magnetic structures, respectively. The Fermi level corresponds to zero energy. 

In order to analyze of the results of the LDA+U calculations in terms of t2 g and e g orbitals one has to chose in the 
local coordinate system (LCS). Usually it is possible to define the LCS where axes are pointed directly along Ti — O 
bonds, but for strongly distorted compounds (like LaTiC>3 and YTiOs) such coordinate system in general might be 
not rectangular. We use the rectangular local coordinate system where the axes are directed as much as possible to 
the nearest oxygens. 

The interorbital on-site Coulomb interaction parameter U and intra-atomic exchange coupling J for ti g shell 
were estimated to be 3.3 eV and 0.8 eV, respectively, using constrained superscell calculations 2 ^ in the framework 
of LMTO method and taking into account screening by e g electrons^ These values are in good agreement with 
previous findings. 29 ! 30 The experimentally observed magnetic structures: YTi0 3 - ferromagnet, LaTi0 3 - G-type 
antifcrromagnet are used for the conventional and constrained LDA+U calculations. 

The top panel of Fig. |2 shows LDA+U electronic structure of LaTiC>3 in details. There are three distinguishable 
sets of bands: completely filled 0-2p bands, partially filled Ti-3d bands and empty La-5o? bands. The bands in the 
energy range from -7.2 eV to -2.7 eV originate mainly from 0-2p states. The gap ^2.3 eV appears between 0-2p and 
the narrow peak of Ti-3c?(t2 g ) states. The band gap in LaTiC>3 is found to be 0.57 eV. It divides Ti-3d(t2g) band into 
two parts in such a way that both the top of the valence and the bottom of the conduction band are predominantly 
formed by Ti-3d(t2g) states (see Fig. 01. The bands lying higher than ^2 eV have basically La-3d, 0-2p and Ti-3d(e 5 ) 
contributions. The calculated local magnetic moment on Ti 3+ ion amounts to 0.78//B being by ~0.2/j,b larger than 
the experimentally observed valued 

The electronic structure obtained in the LDA+U approximation for YTiC>3 is shown in the bottom panel of Fig. [21 
Generally it has similar character. However, the value of the band gap is larger than in LaTi03 (~0.78 eV). The 
value of the local magnetic moment is found to be 0.89/j,s in good agreement with experiment^ 

There are two opposite factors influencing the magnitude of the band gap. On one hand the larger lattice distortions 
in YTi03 due to smaller ionic radii of Y in comparison with La make Ti(3d)-0(2p) hybridization weaker and as a 
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TABLE II: Total energy difference between three magnetic solutions per Ti ion and magnitudes of the band gaps and local 
magnetic moments on Ti ions obtained in LDA+U calculation. Zero energy is the lowest total energy for compound under 
consideration. 
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result the band gap larger. On the contrary, ferromagnetism observed in YTiC>3 leads to an increase of the band 
width resulting in a band gap reduction. 

In order to clarify the role of magnetism we performed additional band structure calculations for magnetic structures 
different from experimental one. 

The FM, AFM-G and AFM-A structures were considered. The values of the band gap and local magnetic moments 
on Ti 3+ ions for all calculated configurations are presented in Tab. |H] One can see that the band gap reduction due 
to ferromagnetism amounts to 21% (0.12 eV) in the case of LaTiOa and 25% (0.26 eV) in YTi03. The modifications 
of an electronic structure with the change of magnetic structure from AFM-G to FM in the vicinity of Fermi level for 
LaTiOs are presented in the inset of FigOU 

Continuing investigation of the interplay between lattice, electronic and magnetic degrees of freedom one can isolate 
the influence of "pure lattice distortions" on the electronic structure of La and Y titanates comparing the results of 
calculations performed in identical magnetic structures. The analysis shows significant changes in the value of the 
band gap. Only due to the lattice distortions it changes by 0.32-0.47 eV depending on the magnetic structure under 
consideration. 

The main reason for such a strong influence of local geometry on the electronic structure of these compounds is 
probably connected with the stronger covalency of d 1 configuration and hence higher sensibility to the distortions in 
comparison with other configurations of rf— shell (with the exclusion of d 9 ). In contrast to the band gap, magnetic 
moment is only reduced by 10% with the change of magnetic configuration. 

In addition, it should be mentioned that our calculations unlike the experiment give for LaTiOs the lowest total 
energy for A-type AFM, although experimentally observed G-type lies quite close. This result is supported by the 
calculation of the exchange interaction parameters^ and connected with the orbital pattern of the compound discussed 
in details in Sec. IVIII 

VI. CALCULATION OF CRYSTAL-FIELD EXCITATION ENERGY IN LDA+U 

The direct study of CFS, for instance by optical measurements, implies the excitation of an electron from one 
energy level to another. Simultaneously with this excitation the external "bath" formed by the rest of the electrons 
can relax giving the system a chance to lower the energy. The simple CFS computations using the centers of gravity 
of corresponding bands, projection or downfolding procedure do not take into account such processes and as a result 
overestimate the value of CFS. In this section we propose the direct calculation of the total energy difference between 
ground and first excited states (CF excitation energy calculation) in ti g shell using the constrain procedure adopted 
for the LDA+U approximation. 



The occupation matrix in LDA+U can be defined in the usual way: 

1 f Ep 

n mm' = / I m G1 nlminlm >{E)dE, (15) 

7T J 

where a is the spin, i denotes the site, n,l,m are principle, orbital and magnetic quantum numbers respectively, 
Ginim inim'(E) — (inlma\(E — H LDA+u )~ 1 \inlm'o-) are the elements of the Green function matrix and H LDA+U is 
a single particle LDA+U Hamiltonian (for its definition see Ref. \2% . The occupation matrix for an isolated ion is 
diagonal. However, for ions in solids it can have more complex structure in low symmetry systems and one needs to 
diagonalize it. The eigenvectors of the i|15[l can be used for transformation of occupation matrix to diagonal form. 
Using this procedure one can obtain the information on the orbitals where the electrons are localized. 
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In fact, in the case of spin-polarized calculations (like LDA+U) there are two occupation matrices for every transition 
metal ion on each site: spin majority and spin minority occupation matrices. According to the Hund's rule, at first 
the spin majority states start to be occupied. Nevertheless, due to the large spatial extension and sizable overlap 
between e g and oxygen p orbitals (see Fig. |3J) there are always some non-negligible occupation numbers for e g states 
in a real band structure calculation for both spins. In order to avoid the influence of these effects (they can lead 
to the "symmetrization" of d-orbitals) we use for the analysis not the occupation matrix for majority spin (Ti is d 1 
system in the compounds under consideration), but the difference between matrices for two spins. Here we suppose 
the oxygen influence on d— states for both spins to be equal. 

The diagonalization of the difference between occupation matrices for majority and minority spins for LaTi03 
gives following eigenvalues: 0, 0, 0.01, 0.02, 0.76. The eigenvector corresponding to the largest eigenvalue defines the 
occupied orbital as a linear combination of all d— orbitals: 

\^ GS )=J2 a i\^)- ( 16 ) 

i 

Thus, in case of LaTiC>3 in the LCS: 

\^% a P° 3 ) = 0.62a;?/ + 0.72yz + 0.32a;z + (17) 
0.02(3y 2 - r 2 ) + 0.02(z 2 - x 2 ). 

Similar decomposition for YTiC>3 is presented in the Table III. 

The excitation energy (from the ground to the first excited state) in this case is the total energy difference between 
states where this orbital is occupied and where it is empty and another one is occupied. According to this scheme 
we performed the calculation where the system is constrained by external potential V cons t r to change the occupied 
orbital: 

V constr = \*GS)SV($GS\- (18) 

In other words V cons tr just pushes up the orbital where the electron has been localized in LDA+U. It is not important 
how big is the correction SV. One needs just be sure that it is large enough to force the electron to hop to another 
orbital. The total energy of excited state does not depend on the value of SV, because the correction is applied to the 
orbital which has to be empty in the excited state. The total energy difference between the ground and first excited 
states is an estimate of the CF excitation energy. 

The results of calculations of CF excitation energies for LaTiC>3 and YTiC>3 are presented in the fourth column 
of Tab. |H They have the same order of magnitude as those obtained using the WF projection in the present work 
and downfolding performed by Pavarini et al& As has been mentioned above, for proper estimation of CF excitation 
energy the relaxation of electron system should be taken into consideration. Such relaxation decreases the energy 
costs of the electron excitation. According to the present results the gain in energy can amount to ~70 meV (~30%). 
However, it is obviously insufficient to reduce the value of CFS significantly. 




FIG. 3: (color online). LaTi03. Ti-3d partial DOS calculated within LDA+U approach. The inset shows Ti-3d partial DOS 
in AFM-G (solid line) and FM (dashed line) configurations. Parts of plots with positive (negative) ordinates denote majority 
(minority) spin DOS. The Fermi level corresponds to zero energy. 
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FIG. 4: (color online) Orbital ordering in LaTiOg (left) and YTi03 (right) in afc-plane obtained in the framework of the LDA+U 
calculations for AFM-G and FM structures respectively. 

Thus, the results of the CFS calculations both using WF projection and constrained LDA+U which takes into 
account the electron system relaxation indicate that the splittings in LaTi03 and YTiC>3 are relatively large. They 
are definitely bigger than the spin-orbital coupling in titanates which is expected to be Aso ~ 20 meV— . The latter 
result is supported by the recent XAS measurements, where shown that the orbital momentum in LaTiC>3 is essentially 
quenched^ 

But the most important conclusion is that CF excitation energy obtained in the present work and value of CFS 
calculated in Ref. as well as extracted from Ti L2.3 XAS^ and optical measurements^ are all the order of 200 meV, 
that makes orbital liquid scenario for LaTiC>3 rather unlikely. 

VII. ORBITAL STRUCTURE 

Orbital degrees of freedom are known to play an important role and should be correctly taken into account in 
theories describing magnetic interactions in strongly correlated materials^* 

As one can see from Tab. IIIII the composition of the occupied orbitals in Y and La titanates is quite different. 
Contribution of xy and yz components to the resulting orbital in LaTi03 are nearly the same and by ~ 40-60 % 
bigger than xz one. There is quite different situation in YTi03, where xz contribution is almost zero. 

This fact can be explained by means of the analysis of the crystal structure distortions, which are quite different in 
these compounds. Distortions in ab— plane in LaTi03 have in general large trigonal D 3 d component. It gives rise to 
square-to-rectangle transformation in this planed and results in localization of the electron on the orbital of almost 
a>ig = {xy + yz + zx) / s/3 character. At the same time the predominant distortions in YTi03 make rhombus from the 
initial square, revealing the sizable contribution of local tetragonal distortions. 



TABLE III: Decomposition of the lowest in energy orbital into xy, xz and yz orbitals. The results are presented in local 
coordinate system where the orthogonal axes point as much as possible to neighboring oxygens. 



LaTiOz YTiOz 



LDA Q.QZxy+Q.myz+OAQxz Q.7Qxy+Qmyz-0.Q6xz 
LDA+U 0.62xy+0.72yz+Q.?,2xz 0.75xy+0.66yz-0.04:Xz 
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However, even in case of LaTi03 the occupied orbital deviates from ai g : 



\{*LDA+u\ai g )\ 2 =91% 
\{^LDA\ai g )\ 2 = 96% 



(19) 
(20) 



The variation of Ti — O distances in YTiOg^ leads to antiferroorbital ordering (Fig. 4,5) inducing according to 
Goodcnough-Kanamori rules the ferromagnetic structure in agreement with experiment 

The basal plane in LaTiC>3 undergoes the elongation in the direction of orthorhombic a— axis. It leads to such kind 
of orbital ordering then the orbitals of Ti ions placed in the ab— plane point almost along the same direction, but not 
to the oxygens (Fig. 01 left panel). Thus, this in-plane "ferroorbital ordering" implies large hopping integrals not only 
between two occupied (AFM interactions), but also between occupied and unoccupied (FM interactions) orbitals 
on different sites. The latter agrees with direct calculation of the exchange interaction parameters in Hartree-Fock 
approximation^ and have to be taken into account in model calculations. 

The presence of a mirror plane in GdFe03-type structure perpendicular to c— axisi£ results in such kind of orbital 
ordering in c— direction that the lobes of the orbitals point to each other (see Fig. [SJ left panel). This is in strong 
contrast to the orbital ordering in ab— plane, where orbitals are directed away from the oxygens. 

The complicated orbital structure of LaTiC>3 does not favor the isotropic magnetic interactions in LaTiC>3. But 
this is not a unique situation. It is generally accepted that YTi03 shows Jahn- Teller distortions and corresponding 
orbital ordering (Ref. and references therein). Experimentally however it also has isotropic magnetic properties, 
which are hard to explain by orbital ordering. This would require an unrealistic set of parameters^ 

Thus the isotropic character of exchange in both LaTiC>3 and YTiC>3 remains an open problem. Whether the orbital 
fluctuations^^ with CF splitting of ~200 meV can resolve this problem, is not clear to us. Another option could 
be that the actual crystal (and consequently orbital) structure of LaTi03 is somewhat different from the accepted 
one - see the recent results in Ref. l38l where the indications of monoclinic distortions have been found. This question 
requires further investigations. 



FIG. 5: (color online) Orbital ordering in LaTiOs (left) and YT1O3 (right) in c-direction obtained in the framework of the 
LDA+ U calculations for AFM-G and FM structures, respectively. 




VIII. CONCLUSIONS 

Following the logic of the present paper we would like to stress two general points. 
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(I) In the framework of LDA+U approximation wc propose the direct method of calculation of crystal-field excitation 
energy. It can be useful for comparison with the results of spectroscopic measurements because it takes into account 
fast relaxations of electronic system. Any electron excitations in spectroscopy are accompanied by such relaxations 
and it should be definitely taken into consideration in the calculation of these excitations. 

Moreover, this scheme is more reliable then others because it uses not eigenvalues of LDA but its total energies. In 
addition, the method does not have an ambiguity in the definition of the basis set as any projection or downfolding 
procedure does, because all states included in the Hamiltonian arc used. The latter is especially important for low 
symmetry systems where it is not clear what kind of linear combination of d— wave function and in which coordinate 
system should be taken for construction of a few orbital tight-binding Hamiltonian. 

(II) We probe this method on Y and La titanates, for which the value of crystal- field excitation energy is indeed 
essential. Using Wannier function projection of LDA full-orbital Hamiltonian it is found that CFS between the lowest 
and the next in energy t2 g level is 230 meV and 180 meV for La and Y titanates, in agreement with previous estimated 

According to the results of present calculations, fast electronic relaxations reduce CFS by ~30% (it amounts to 
70 meV in case of LaTiOs) as compared with the difference of LDA one electron energies or CFS obtained by 
Madelung point charge summation^ However, the value of crystal-field excitation energy is still large enough to make 
orbital liquid formation in titanates rather unlikely. At the same time the intricate orbital ordering in contrast to the 
expectations from the model calculations does not favor isotropic magnetic interaction in LaTiOs. 
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